restart;
with(LinearAlgebra):


#3D steady velocity field U3 in first-order Taylor approximation
uuf := uu + uu_x*x + uu_y*y + uu_z*z;
vvf := vv + vv_x*x + vv_y*y + vv_z*z;
wwf := ww + ww_x*x + ww_y*y + ww_z*z;
U3 := Vector([uuf,vvf,wwf]);


#current particle velocity
V3 := Vector([u,v,w]);

#gravity Vector
G3 := Vector([ug,vg,wg]);


#3D Spatial Jacobian matrix J3:
J3 := Matrix([
[diff(uuf,x),diff(uuf,y),diff(uuf,z)],
[diff(vvf,x),diff(vvf,y),diff(vvf,z)],
[diff(vvf,x),diff(wwf,y),diff(wwf,z)]]);


#6D vector field U6
dx := u;
dy := v;
dz := w;
du := (uuf-u)/r + ug;
dv := (vvf-v)/r + vg;
dw := (wwf-w)/r + wg;
U6 := Vector([dx,dy,dz,du,dv,dw]);


#6D Jacobian matrix J6
J6 := Matrix([
[diff(dx,x),diff(dx,y),diff(dx,z),diff(dx,u),diff(dx,v),diff(dx,w)],
[diff(dy,x),diff(dy,y),diff(dy,z),diff(dy,u),diff(dy,v),diff(dy,w)],
[diff(dz,x),diff(dz,y),diff(dz,z),diff(dz,u),diff(dz,v),diff(dz,w)],
[diff(du,x),diff(du,y),diff(du,z),diff(du,u),diff(du,v),diff(du,w)],
[diff(dv,x),diff(dv,y),diff(dv,z),diff(dv,u),diff(dv,v),diff(dv,w)],
[diff(dw,x),diff(dw,y),diff(dw,z),diff(dw,u),diff(dw,v),diff(dw,w)]]);

x := 0;
y := 0;
z := 0;

#compare (9) in the paper:
J6;

#compare (11) in the paper:
W := (U3-V3)/r + G3;

#compare (12) in the paper:
Multiply(J6,U6);